!
!     TEST OF THE SUBROUTINE LUNEPH.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,m,o-z)
      dimension r(3)
!
! --- Parameters of computation ----------------------------------------
!
!     t0 : Initial Julian date.
!     dt : Time interval (day).
!     nd : Number of dates.
!      
      parameter (t0=2440222.5d0,dt=2000.0d0,nd=8)
!           
! --- Files ------------------------------------------------------------
!
!     11 : Polynomials :  luneph.dat
!     12 : Results     :  luneph.txt
!      
      open (11,file='luneph.dat',status='old')
      open (12,file='luneph.txt')
!
! --- Computation ------------------------------------------------------
!
      write (12,1001) t0,dt,nd      
      do n=1,nd
         tdj=t0+(n-1)*dt
         call LUNEPH (tdj,11,r,ierr)
         if (ierr.ne.0) then
            write (*,"(/2x,'ERREUR LUNEPH : ',i2/)") ierr
            stop
         endif
         write (12,1002) n,tdj,r
         write (*,1002)  n,tdj,r
      enddo
!
      stop
!
! --- Formats ----------------------------------------------------------
!
1001  format (/2x,'LUNAR EPHEMERIS : LUNEPH'/
     .        /2x,'INITIAL JULIAN DATE  : ',f9.1
     .        /2x,'TIME INTERVAL (DAY)  : ',f9.4
     .        /2x,'NUMBER OF DATES      : ',i9/)
1002  format (2x,i4.4,f15.5,
     .        2x,f12.6,' deg',2x,f10.6,' deg',2x,f14.6,' km')
!
      end
!
!
!
      subroutine LUNEPH (tdj,nul,r,ierr)
!-----------------------------------------------------------------------
!
!     Ref : 2003/01                     
!
! --- Object -----------------------------------------------------------
!
!     Lunar Ephemeris LUNEPH. 
!     Geocentric coordinates : Longitude, Latitude, Distance. 
!     Reference frame : dynamical ecliptic and equinox of J2000.0.
!
! --- Input ------------------------------------------------------------
!
!     tdj         Julian date (real*8).
!
!     nul         Logical unit number of LUNEPH file (integer).
!
! --- Output -----------------------------------------------------------
!
!     r(3)        Coordinates table (real*8).
!                 r(1) : Longitude (degree, between 0 and 360).
!                 r(2) : Latitude  (degree).
!                 r(3) : Distance  (km).
!
!     ierr        Error index (integer).
!                 0 : No error.
!                 1 : Reading file error. 
!                 2 : Date error.
!
! --- Remark -----------------------------------------------------------
!
!     The LUNEPH polynomials file 'luneph.dat' has to be opened before 
!     the first call to LUNEPH subroutine.
!
! --- Parameters of Chebychev polynomials ------------------------------
!
!     t0 : Initial Julian date in the period 1950-2050.
!     dt : Lengh of time interval of validity of the polynomials (day).
!     ni : Number of time intervals in the period.
!     nc : Number of coefficients per coordinate.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      logical first
!
      dimension r(3)
      dimension cf(9223,12,3),tn(12)
!      
      data first/.true./
!
      save 
!
! --- Reading the coefficients of the polynomials LUNEPH ---------------
!
      if (first) then
         write (*,'(/a/)') ' << READING THE COEFFICIENTS LUNEPH >> '
         read (nul,1000,iostat=nerr) t0,dt,ni,nc
         do n=1,ni
            read (nul,1001,iostat=nerr)
            if (nerr.ne.0) then
               ierr=1
               return
            endif   
            do j=1,nc
               read (nul,1002,iostat=nerr) (cf(n,j,k),k=1,3)
               if (nerr.ne.0) then
                  ierr=1
                  return
               endif   
            enddo
         enddo
         first=.false.
         tf=t0+ni*dt
      endif
!
! --- Check of validity of dates ---------------------------------------
!
      if (tdj.lt.t0.or.tdj.gt.tf) then
         ierr=2
         return
      else     
         ierr=0
      endif   
!
! --- Substitution of time in the polynomials --------------------------
!
      n=(tdj-t0)/dt+1
      if (n.gt.ni) n=ni
      ti=t0+(n-1)*dt
      x=2.d0*(tdj-ti)/dt-1.d0
      tn(1)=1.d0
      tn(2)=x
      w=x+x
      do i=3,nc
         tn(i)=w*tn(i-1)-tn(i-2)
      enddo
      do j=1,3
         r(j)=0.d0
         do i=1,nc
            r(j)=r(j)+cf(n,i,j)*tn(i)
         enddo
      enddo
!
! --- Reduction of the longitude between 0 and 360 degrees -------------
!
      r(1)=mod(r(1),360.d0)
      if (r(1).lt.0.d0) r(1)=r(1)+360.d0
!
      return
!
! --- Formats ----------------------------------------------------------
!
1000  format (1x,f9.1,2x,f4.1,2x,i5,2x,i2)      
1001  format (1x)
1002  format (3f20.12)
!
      end
